#ifndef CSB_BELLE_HPP_
#define CSB_BELLE_HPP_

#include <iostream>

#include <csb-util.hpp>

#include "particle/Particle.h"
#include "particle/PID.h"
#include "particle/utility.h"
#include "event/BelleEvent.h"
#include "basf/module.h"
#include "basf/module_descr.h"
#include "benergy/BeamEnergy.h"
#include "kid/atc_pid.h"
#include "ip/IpProfile.h"
#include "helix/Helix.h"
#include "kfitter/kvertexfitter.h"

#include <mdst/mdst.h>

#include "panther/panther.h"
#include MDST_H
#include MDST_OBS_H
#include EVTCLS_H
#include BELLETDF_H
#include RUN_INFO_H

#define LGNDR2(X) ((3.*X*X-1.)/2.)
#define LGNDR3(X) ((5.*X*X*X-3.*X)/2.)
#define LGNDR4(X) ((35.*X*X*X*X-30.*X*X+3.)/8.)
#define LGNDR5(X) ((63.*X*X*X*X*X-70.*X*X*X+15.*X)/8.)
#define LGNDR6(X) ((231.*X*X*X*X*X*X-315.*X*X*X*X+105.*X*X-5.)/16.)
#define LGNDR7(X) ((429.*X*X*X*X*X*X*X-693.*X*X*X*X*X+315.*X*X*X-35.*X)/16.)
#define LGNDR8(X) ((6435.*X*X*X*X*X*X*X*X-12012.*X*X*X*X*X*X+6930.*X*X*X*X-1260.*X*X+35)/128.)

static const double E_MASS = 0.000511;
static const double CHARGED_MASS[5] = {0.000511, 0.105658366, 0.13957018,
				       0.493667, 0.939272046};


#ifdef BELLE_NAMESPACE
namespace Belle {
#endif //BELLE_NAMESPACE

  void track_combinations(std::vector<std::vector<Mdst_charged> >& output, 
			  std::vector<Mdst_charged> pos, 
			  std::vector<Mdst_charged> neg, int npos, int nneg) 
  {
    if (npos > pos.size())
      return;

    if (nneg > neg.size())
      return;

    std::vector<std::vector<Mdst_charged> > pos_combos; 
    csb::unordered_combinations<Mdst_charged>(pos, pos_combos, npos);

    std::vector<std::vector<Mdst_charged> > neg_combos;
    csb::unordered_combinations<Mdst_charged>(neg, neg_combos, nneg);

    for (std::vector<std::vector<Mdst_charged> >::iterator ipos_combos=pos_combos.begin(); 
	 ipos_combos != pos_combos.end(); ipos_combos++) 
      {	
	for (std::vector<std::vector<Mdst_charged> >::iterator ineg_combos=neg_combos.begin();
	     ineg_combos != neg_combos.end(); ineg_combos++)
	  {
	    std::vector<Mdst_charged> tmp;
	    if (npos >= nneg) {
	      for (std::vector<Mdst_charged>::iterator ipos = ipos_combos->begin(); 
		   ipos != ipos_combos->end(); ipos++)
		tmp.push_back(*ipos);
	      
	      for (std::vector<Mdst_charged>::iterator ineg = ineg_combos->begin(); 
		   ineg != ineg_combos->end(); ineg++)
		tmp.push_back(*ineg);
	    }
	    else {	      
	      for (std::vector<Mdst_charged>::iterator ineg = ineg_combos->begin(); 
		   ineg != ineg_combos->end(); ineg++)
		tmp.push_back(*ineg);

	      for (std::vector<Mdst_charged>::iterator ipos = ipos_combos->begin(); 
		   ipos != ipos_combos->end(); ipos++)
		tmp.push_back(*ipos);

	    }	      

	    output.push_back(tmp);
	  }
      }

  }

	      
		



  long two_body_idhash(long p1id, long p2id) {
    long ids[2] = {p2id, p1id};
    
    long h = p2id;
    h += 1000 * p1id;
    
    return h;
  }

  long three_body_idhash(long p1id, long p2id, 
			 long p3id) 
  {
    long ids[3] = {p3id, p2id, p1id};
  
    long h = 0;
    
    long factor = 1000;
    for (int i=0; i < 3; i++) 
      h += ids[i] * (long) pow(factor, i);    
    
    return h;
  }

long four_body_idhash(long p1id, long p2id, 
		      long p3id, long p4id) {
  long ids[4] = {p4id, p3id, p2id, p1id};
  
  long h = 0;
  
  long factor = 1000;
  for (int i=0; i < 4; i++) 
    h += ids[i] * (long) pow(factor, i);    
  
  return h;
}

HepLorentzVector initial_state() {
    double E_BEAM = BeamEnergy::E_beam_corr();
    double E_HER = BeamEnergy::E_HER();
    double E_LER = BeamEnergy::E_LER();
    double cross_angle = BeamEnergy::Cross_angle();

    HepLorentzVector p4_her(0.0,
			    0.0,
			    +sqrt(E_HER*E_HER-E_MASS*E_MASS),
			    E_HER);
    p4_her.rotateY(cross_angle);
    HepLorentzVector p4_ler(0.0,
			    0.0,
			    -sqrt(E_LER*E_LER-E_MASS*E_MASS),
			    E_LER);
    
    HepLorentzVector initial = p4_ler + p4_her;

    return initial;
}
  

  HepLorentzVector charged_track_p4(const HepPoint3D& IP, 
				    const Mdst_charged& trk, 
				    int type) 
  {
    Mdst_trk_fit& trkFit = trk.trk().mhyp(type); 
	  
    HepVector a(5,0);
    a[0] = trkFit.helix(0);
    a[1] = trkFit.helix(1);
    a[2] = trkFit.helix(2);
    a[3] = trkFit.helix(3);
    a[4] = trkFit.helix(4);
    
    HepPoint3D pivot(trkFit.pivot(0), trkFit.pivot(1), trkFit.pivot(2));
    Helix ltrk(pivot, a);
    //	  ltrk.pivot(HepPoint3D(0.,0.,0.));
    ltrk.pivot(IP);
    
    return  ltrk.momentum(0., CHARGED_MASS[type]);
  }
   
  
  
 
/*
 .......................................................................
 . FWJET2 - Subroutine to compute event sphericity 
 . 
 . Inputs    : NTRK -- number of tracks in momentum array 
 .           : PTRK -- array dimensioned (3,NTRK) containing 3-momenta 
 .           :         of tracks in the event 
 .           : EBEAM -- beam energy of the event 
 . Outputs   : H0 -- 0th moment of the event 
 .           : R1 -- H1/H0 of the event 
 .           : R2 -- H2/H0 of the event 
 .           : R3 -- H3/H0 of the event 
 .           : R4 -- H4/H0 of the event 
 . 
 . Calls     : LEGN(X) -- INTERNAL functions for Legendre polynomials 
 . Called    : Anyone who wants the Fox/Wolfram moments of the event 
 . Author    : F. Morrow  09/02/89  14.47.58 
 .             Transfer to C by kfjack 2/9,1999
 . 
 . Detailed description 
 .     FWJET2 determines the Fox/Wolfram moments of the event, up to 
 .     order 5. 
 .......................................................................
*/

void fwjet2(int NTrk, float *PTrk, float EBeam, 
	    float *H0, float *R1, float *R2, float *R3, float *R4)
{
  int i, j;
    float sss ,cxang, ppari, pparj, delh0;
    float h, h1, h2, h3, h4;
    
    if (NTrk >= 2) {
      
      if (EBeam < 4.0 || EBeam > 6.5) EBeam = 5.29;
      sss = EBeam*EBeam*4;
    
      h=h1=h2=h3=h4=0;
      
      for (i=0;i<NTrk;i++) {
	ppari = sqrt(PTrk[i*3]*PTrk[i*3]+
		     PTrk[i*3+1]*PTrk[i*3+1]+
		     PTrk[i*3+2]*PTrk[i*3+2]);
	
	for (j=i;j<NTrk;j++) {
	  
	  if (i!=j) {
	    pparj = sqrt(PTrk[j*3]*PTrk[j*3]+
			 PTrk[j*3+1]*PTrk[j*3+1]+
			 PTrk[j*3+2]*PTrk[j*3+2]);
	    
	    cxang = (PTrk[i*3]*PTrk[j*3]+ 
		     PTrk[i*3+1]*PTrk[j*3+1]+
		     PTrk[i*3+2]*PTrk[j*3+2])/(ppari * pparj);
	    
	    delh0 = 2*ppari*pparj;
	    
	    h += delh0;
	    h1 += delh0 * cxang;
	    h2 += delh0 * LGNDR2(cxang);
	    h3 += delh0 * LGNDR3(cxang);
	    h4 += delh0 * LGNDR4(cxang);
	  } else {
	    delh0 = ppari*ppari;
	    h += delh0;
	    h1 += delh0;
	    h2 += delh0;
	    h3 += delh0;
	    h4 += delh0;
	  }
	}
      }
    }
    
    if (h > 1e-7f) {
      *H0 = h / sss;
      *R1 = h1 / h;
      *R2 = h2 / h;
      *R3 = h3 / h;
      *R4 = h4 / h;
    }else {
      *H0=*R1=*R2=*R3=*R4=-1.0;
    }
  }

  float calc_R2() {
    Mdst_charged_Manager& m_charged = Mdst_charged_Manager::get_manager();
    Mdst_gamma_Manager& m_gamma = Mdst_gamma_Manager::get_manager();
        
    int n_trk = m_charged.count() + m_gamma.count();
    float H0, R1, R2, R3, R4;

    float p_trk[n_trk*3];

    float beam_e = BeamEnergy::E_beam_corr();

    int trkno = 0;
    for (Mdst_charged_Manager::iterator it = m_charged.begin(); 
	 it != m_charged.end(); it++, trkno++)
      {
	p_trk[3*trkno]     = it->px();
	p_trk[3*trkno + 1] = it->py();
	p_trk[3*trkno + 2] = it->pz();
      }

    for (Mdst_gamma_Manager::iterator it = m_gamma.begin();
	 it != m_gamma.end(); it++, trkno++)
      {
	p_trk[3*trkno]     = it->px();
	p_trk[3*trkno + 1] = it->py();
	p_trk[3*trkno + 2] = it->pz();
      }	

    fwjet2(n_trk, p_trk, beam_e, &H0, &R1, &R2, &R3, &R4);

    return R2; 
  }



#ifdef BELLE_NAMESPACE
}
#endif


#endif // CSB_BELLE_HPP_
